Finite Larmor radius effects on non-diffusive tracer transport in a zonal flow 
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Finite Larmor radius (FLR) effects on non-diffusive transport in a prototypical zonal flow with 
drift waves are studied in the context of a simplified chaotic transport model. The model consists 
of a superposition of drift waves of the linearized Hasegawa-Mima equation and a zonal shear flow 
perpendicular to the density gradient. High frequency FLR effects are incorporated by gyroaveraging 
the E X B velocity. Transport in the direction of the density gradient is negligible and we therefore 
focus on transport parallel to the zonal flows. A prescribed asymmetry produces strongly asymmetric 
non- Gaussian PDFs of particle displacements, with Levy flights in one direction but not the other. 
For k±pth ~ 0, where k± is the characteristic wavelength of the flow and pth is the thermal Larmor 
radius, a transition is observed in the scaling of the second moment of particle displacements, 
a^ ^ t"' . The transition separates ballistic motion, 7 ~ 2, at intermediate times from super- 
diffusion, 7 = 1.6, at larger times. This change of scaling is accompanied by the transition of the 
PDF of particle displacements from algebraic decay to exponential decay. However, FLR effects seem 
to eliminate this transition. In all cases, the Lagrangian velocity autocorrelation function exhibits 
non-diffusive algebraic decay, C ~ t~^, with (^ = 2 — 7 to a good approximation. The PDFs of 
trapping and flight events show clear evidence of algebraic scaling with decay exponents depending 
on the value of k±pth- The shape and spatio-temporal self-similar anomalous scaling of the PDFs 
of particle displacements are reproduced accurately with a neutral, a — P, asymmetric effective 
fractional diffusion model where a and /3 are the orders of the spatial and temporal fractional 
derivatives. 

PACS numbers: 52.25.Gj,52.35.Kt,52.65.Cc,05.40.Fb,05.45.Pq,52.25.Fi,52.65.-y 



I. INTRODUCTION 

Plasma turbulence presents a challenge to multiscale 
models of transport in applications such as magnetic fu- 
sion confinement, stellar accretion disks and galactic dy- 
namos. Simulations of turbulent transport involve non- 
linear interactions at disparate scales, which often makes 
numerical computations expensive and analytic methods 
intractable. As an alternative, one may consider models 
of intermediate complexity that incorporate important 
aspects of transport within a relatively simple reduced 
description. In this paper we follow this approach and 
present a numerical study of the role of finite Larmor ra- 
dius (FLR) effects on non-diffusive poloidal transport in 
zonal shear flows using a reduced E x B Hamiltonian test 
particle transport model. 

Following Ref. [l|, we model the flow as a superposi- 
tion of a shear flow and drift waves obtained from the lin- 
earized Hascgawa-Mima (HM) equation 0] ■ Test particle 
characteristics in this flow are generally not integrable 
and exhibit chaotic advection, also known as Lagrangian 
turbulence, which reproduces key ingredients of particle 
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transport in more complex flows. High frequency FLR 
effects are incorporated by solving the test particle equa- 
tions of motion for the gvroaveraged E x B velocity. As 
demonstrated by Ref. [3|, we compute the gyroaverage 
using a discrete A'^-polygon approximation. 

We adopt a statistical approach and apply non- 
diffusive transport diagnostics to large ensembles of par- 
ticles. One of the simplest diagnostics is the scaling of 
the second moment of particle displacements, cr^(t) = 
{[Sy — {5y)Y), where 6y = 6y{t) denotes the particle's 
displacement and ( ) denotes the ensemble average. In 
the standard diffusion case, cr^(t) ^ t, linear scaling al- 
lows the definition of an effective diffusivity as the ra- 
tio Deff — CT^(i)/(2i) in the limit of large t. However, 
in the case of non-diffusive transport, cr'^(t) ~ f^ with 
7 7^ 1. When < 7 < 1, the growth of the variance 
is slower than diffusion and transport is sub-diffusive. 
When 1 < 7 < 2 transport is super- diffusive, which 
means the spreading is faster than diffusion, and the dis- 
placements may be Levy flights !4|. In both super- and 
sub-diffusion, characterization of transport as a diffusive 
process with an "effective diffusivity" -De// breaks down 
because D^fj — > when < 7 < 1, and -De// — ^ 00 
when 1 < 7 < 2. Other measures of non-diffusive trans- 
port, which will be discussed in detail later, include non- 
Gaussianity of the probability distribution of displace- 
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FIG. 1: (Color online) Contour plots of electrostatic potential 
(j>. Panel (a) shows a snapshot of the potential obtained from 
a direct numerical simulation of the Hasegawa-Mima equa- 
tion ([5]). Panel (b) shows (^ at a fixed time according to the 
chaotic Hamiltonian transport model in Eq. ^. The thick 
line limiting the central vortex in (b) is the separatrix. Par- 
ticles inside the separatrix are trapped, and, as the arrows 
show, particles outside the separatrix are transported by the 
zonal flow. The Hamiltonian model in (b) provides a reduced 
description of E x B transport dominated by vortices and 
zonal flows as highlighted by the rectangle in (a). 



ments (propagator), slow decay of the Lagrangian veloc- 
ity autocorrelation function, the presence of long jumps 
(Levy flights) and long waiting times, and the non-local 
{i.e., non-Fickian) dependence of fluxes on gradients. A 
general review of non-diffusive transport can be found in 
Ref. [^], and discussions focusing on plasmas can be 
found in Ref. [IH. 

Test particle transport in HM flows, as in Fig.[lja), has 
been studied in Refs. P, i, i, M, EH, iB, [11 ■ In Ref. Q , 
which did not include FLR effects, it was shown that 
zonal flows give rise to Levy flights and strongly asym- 
metric non-Gaussian PDFs of particle displacements. 
References [3, [l0| addressed the role of FLR effects but 
restricted attention to diffusive transport. More recently, 
Ref. [13[ considered FLR effects in non-diffusive trans- 
port in HM turbulence and concluded that the expo- 
nent 7 does not change appreciably with the Larmor ra- 
dius but that the effective diffusion coefficient is reduced. 
There is a very close connection between drift waves as 
described by the HM equation and Rossby waves as de- 
scribed by the quasigeostrophic equation, see for exam- 
ple Ref. [lj|. Therefore, statistical test particle studies 
in fluid mechanics, such as Refs. |15l.ll6j. are in principle 
applicable to drift wave transport. 

The main new results presented here, which to our 
knowledge have not been reported in the literature be- 



fore, include: (i) a transition from algebraic to exponen- 
tial decay in the tails of PDFs of particle displacements 
accompanied by a transition from ballistic (7 « 2) to 
super- diffusive (1 < 7 < 2) transport; (ii) a numeri- 
cal study of the role of FLR on the Lagrangian veloc- 
ity autocorrelation function and on the particle trap- 
ping and particle flight PDFs; (iii) the construction of 
a effective fractional diffusion model that reproduces the 
shape and the spatio-temporal anomalous self-similar 
scaling of the PDF of particle displacements. In re- 
cent years, fractional diffusion models have been ap- 
plied to describe non-diffusive plasma transport, e.g. 
Refs. [13, M, M, M,M,^- Although the present work 
focuses on a prototypical model of transport, the diagnos- 
tics used and the non-diffusive phenomenology discussed 
here might be of relevance to the study of transport in 
more general flows dominated by coherent structures like 
zonal flows and eddies. Despite the fact that these co- 
herent structures are ubiquitous in simulations and ex- 
periments [ijj [23, [231, their influence on non-diffusive 
transport is not well understood. In this regard, Ref. |25j 
showed evidence of non-diffusive transport in gyrokinetic 
turbulence for "intermediate" simulation times. 

The rest of the paper is organized as follows. In 
Sec. [H] the E x B transport model with and without 
FLR effects is explained. Section Hill shows a benchmark 
of the numerical method against an exact solution for 
the particle propagator in a parallel flow. Section IIVI 
presents a summary of Lagrangian diagnostics to study 
non-diffusive transport. The main numerical results are 
presented in Sec. El Section IVj describes the anoma- 
lous self-similarity properties of the PDF of particle dis- 
placements and presents an effective fractional diffusion 
model. Section FVlII contains the conclusions. 



II. TRANSPORT MODEL 

We follow a Lagrangian approach to study transport 
and consider large ensembles of discrete particles moving 
in a prescribed flow. We limit attention to test particles, 
neglecting self-consistency effects and assuming that the 
particles are transported by the flow without modifying 
it. When finite Larmor radius (FLR) effects can also be 
neglected, the dynamics are determined by a drift equa- 
tion which, in the E x B approximation, is 

dr _ E X B 

where r = {x, y) denotes the particle position, E is the 
electrostatic field, and B is the magnetic field. Writing 
B = Bqz, and E = — V0(a;, y, t), Eq. (P) can be equiva- 
lently written as the Hamiltonian dynamical system 

dx d(j) dy d(j) 

dt dy ' dt dx ' 

where the electrostatic potential is analogous to the 
Hamiltonian, and the spatial coordinates are the canon- 
ical conjugate phase space variables. 



(2) 



For relatively high energy particles or for a flow vary- 
ing relatively rapidly in space, the zero Larmor radius ap- 
proximation fails and it is necessary to incorporate FLR 
effects. A simple, natural way of doing this is to substi- 
tute the E X B flow on the right hand side of Eq. ([2]), 
which is evaluated at the location of the guiding center, 
by its value averaged over a ring of radius p, where p is 
the Larmor radius [3[. Formally, the procedure is given 

by 



takes the form 



dx 



90 
dy 



dy 
dt 



dx 



(3) 



where the gyroaverage, O^, is defined as 



1 r^- 

(*)e = — / * (a- + /ocos6', j/-|-psin6')(i6'. (4) 

This is a good approximation provided the gyrofrcqucncy 
is greater than other frequencies in the system. 

In the HM model for drift waves the electrostatic po- 
tential is determined from [3| 



[dt + (z X V0) • V] (V^ 



/3a;) = , (5) 



where the x coordinate corresponds to the direction of the 
density gradient driving the drift-wave instability, and y 
corresponds to the direction of propagation of the drift- 
waves. In toroidal geometry, x is analogous to a nor- 
malized coordinate along the minor radius, and y is a 
poloidal-like coordinate. Here we assume a slab approx- 
imation and treat {x, y) as Cartesian coordinates. The 
parameter /? = no(x)'/no(x) measures the scale length 
of the density gradient. We model the electrostatic po- 
tential (test particle Hamiltonian) as a superposition of 
an equilibrium zonal shear flow, ipo^x), and the corre- 
sponding eigenmodes of Eq. (O, (pj{x), with perpendic- 
ular wave numbers, kj_j, and frequencies, Cjk±j, 

N 

(j) = (Pq{x) +^ Ej ipj{x) COS k±j{y - Cjt) . (6) 
i=i 

We consider a monotonic zonal flow of the form 



''^yfiix) = tanh(a;) . 



(7) 



In this case, depending on the parameter values, there is 
a band of unstable modes bounded by two regular neutral 
modes with eigenfunctions []| 



ifj = [1 + tanhx] 2 [1 — tanhx] 



(8) 



Since these modes arc neutral, ci and C2 are real and 
the corresponding values of kj_j are obtained from the 
linear dispersion relation. Neutral modes are important 
because they describe dynamics near marginal stability. 
Following Ref. [l|, we consider a traveling wave pertur- 
bation of the first neutral mode. The electrostatic poten- 
tial in the co-moving reference frame of the neutral mode 



(f) = hi (coshx) + ipi{x) [ei cos k±iy + 
£2 coa{k^2y - ^t)] - CiX. 



(9) 



The first term on the right hand side of Eq. ^ is the po- 
tential of the shear flow in Eq. ([7]), and lo is the frequency 
of the perturbation. The wavcnumbers perpendicular to 
the uniform magnetic field, k±i and kj_2, characterize 
the size of E x B eddies, while £i and £2 give the ampli- 
tudes of the waves. When computing k±pth to compare 
the scale length of the eddies in this flow to the thermal 
gyroradius, we use the mean value kj_ ~ {kj_i + k±2)/2. 
When £2 = the Hamiltonian in Eq. ([9|) is time in- 
dependent, and the test particles follow contours of con- 
stant (/) shown in Fig. [Ijb). In this case, particles inside 
the separatrix remain trapped and those outside the sepa- 
ratrix are always untrapped with y > left of the vortices 
and y < right of the vortices. However, when there 
is a time dependent perturbation, i.e. when £2 7^ in 
Eq. ([9]), the E x B particle trajectories are in general not 
integrable. In this case, the separatrix breaks and forms a 
stochastic layer where test particles alternate chaotically 
between being untrapped in the zonal flow and being 
trapped inside the vortices. This is the phenomenon of 
chaotic transport that has been studied in both plasmas 
and fluid systems, see for example Refs. [1, [l^, [2^, \2l\ 
and references therein. As Fig. [Ija) illustrates, the sim- 
ple Hamiltonian model in Eq. Q provides a reduced de- 
scription of E X B eddies embedded in a background zonal 
flow in HM turbulence. 



III. NUMERICAL METHOD 

The zero Larmor radius calculations are based on the 
Hamitonian-like equations of Eq. ^. For the numerical 
integration of these equations we used the second-order 
symplectic predictor-corrector scheme of Ref. [28[ with a 
fixed time step of 0.05 and 8 iterations in the predictor- 
corrector loop. These parameters were chosen based on 
numerical convergence studies and by monitoring the ac- 
curacy of energy conservation. For the model parame- 
ters we used £1 = 0.5, £2 = 0.2, ci = 0.4, k±i = 6.0, 
k±2 = 5.0 and uj = 6.0. This choice is motivated by 
Refs. [il, [I^ where it was shown that, for this set of 
parameters, test particles exhibit strongly asymmetric, 
non-Gaussian statistics. As such, these parameters are 
a good starting point to study the role of FLR effects 
on non-diffusive transport. For the initial conditions 
we used an ensemble of particles located in the vicin- 
ity of the hyperbolic fixed point of the Hamiltonian at 
{xo,yo) ~ (—1,-0.5). This localization guarantees that 
a large fraction of the particles will stay in the stochastic 
layer and undergo chaotic transport. Other choices of 
initial positions can lead to integrable motion with par- 
ticles permanently either inside the eddies, circling, or 
outside, following the zonal flow. 



The only difference between the zero and finite Larmor 
radius calculations is in the evaluation of the velocity of 
the test particle. Assuming fast gyration in a strong B 
field, the gyroaverage of the E x B velocity is computed 
over a circle of radius p, where p is the Larmor radius 
of the particle. Throughout this paper we will assume a 
Maxwellian equilibrium distribution for the Larmor radii 
of the test particles of the form 



H{p) = -^e 

Pth 



-p'ipiu 



(10) 



normalized according to /„ H{p)pdp — 1. For the nu- 
merical computation of the gyroaverage we approximate 
the circle with an inscribed polygon with A^g-sides and 
approximate the integral over the circle as the average 
over the vertices of the polygon. This method, widely 
used in kinetic particle codes {e.g. [3|), simply samples 
the field on the gyration arc at a small number of equally 
spaced points. For example, the 8-point (octagon) ap- 
proximation evaluates the gyroaverage by considering 
Ng ~ 8 points distributed around the circle in equal in- 
crements, i.e., at 9 = {2tt/8,2tt/7, . . .2tt}. If the mean 
gyroradius, (p) — {y/n/2)pth, becomes large relative to 
the typical scale length, ~ l/fc±, of the flow, i.e., if 
k±pth ^ 1, the number of points used to compute the 
gyroaverage must be increased to maintain the same level 
of accuracy. 

The error involved in the approximation of the gyroav- 
erage on Ng for a given value of k±pth and, therefore, a 
benchmark for the accuracy of the numerical scheme can 
be studied by considering the following parallel flow in 
arbitrary geometry 



(f> — (f>o cos{k±x) . 



(11) 



The main object of interest is the probability distribution 
function of particle displacements, or propagator, P = 
P{y, t\y' , t'), which gives the probability for a particle to 
be at y' at time t' if it was at y at time t. Since Vx = Q for 
this choice of (j), we restrict study to the y direction. The 
function P depends on kj_pth and the goal is to study 
the error in the numerical evaluation of P as function of 
k±pth and the value of Ng used in the approximation of 
the gyroaverage. As discussed in Appendix [XI the exact 
propagator for Eq. pTJ) is given by 



P{y,t\y',t') 



1 



Uoit-t') 



77^(0, C 



1 {y~y') 



Uo [t - 1') 



with 



^(0 = 



■E 



-{z,/k^pth? 



ik±Pthf tt \MZ'^)\ 



(12) 



(13) 



where Zi = Zi(Q denotes the i-th zero of the equation 
Jo{zi) — C = 0- Here, Jq is the order zero Bessel function 
of the first kind. For a given (, the number of zeros of 
this equation is Nz which goes to oo as ^ goes to zero. 




FIG. 2: (Color online) Particle propagator for finite Lar- 
mor radius transport in the parallel shear flow of Eq. (|ll|l . 
Panel (a) corresponds to k±pth = 3.0 and (b) corresponds to 
k±pth ~ 5.0. The solid line denotes the exact analytical result 
in Eq. (|12p . the dashed line and the marked line (shown only 
in (b)) denote the 8-point and the 16-point average numerical 
results, respectively. 



Note also that because the minimum and maximum val- 
ues of Jo are —0.4025 and 1, respectively, no zero ex- 
ists for C < -0.4025 or ^ > 1- Therefore, P identically 
vanishes outside the interval C € (—0.4025,1). Despite 
its apparent complexity, this analytical result provides a 
valuable benchmark to assess the accuracy of the gyroav- 
erage computation. 

Figure [2] compares the exact propagator in Eq. ([72]) 
with the propagator obtained from direct numerical inte- 
gration of the gyroaverage equations of motion in Eq. ([3|) 
for different values of k±pth and Ng. The FLR effects 
significantly change the k±pth ~ propagator, which is 
a (5-function centered at C = 1: PiyiAv' ^^') = (l/f^o(* — 
t'))5[C - 1). It is observed that for ki_pth = 3.0, Ng^% 
produces relatively good results, although it misses the 
small spike in P around Sy/Uot ^ 0.25. Other Ng = 8 
cases with k±pth < 3.0 (not shown) give nearly exact 
agreement. However, for k±pth = 5.0, the Ng = 8 av- 
erage departs significantly from the exact result. This 
failure means that choosing Ng > 8, such as Ng = 16, 
is necessary. One is led to conclude that the Ng = 8 
method should not be used for values of k±pth ^ 3.0. 

This statement is further supported by an assessment 
of accuracy when representing Jo((-) as a finite sum based 
on the integral 



1 /-^^ 



cos(t sin T)dT . 



(14) 



The Bessel function is used in spectral simulations of the 
gyrokinetic equation, which gives the spectral technique 
an advantage that we cannot use here. The Bessel in- 
tegral representation may be discretized and evaluated 



using different numbers of terms in tfie sum. Additional 
terms in the sum reduce the error of discretization just as 
increasing Ng reduces the error of discrete gyroaveraging. 
When the integral is approximated with 8 or 16 equally 
spaced points between and 2tt, the result agrees to 0.1% 
with the value of Jaii) up to i = 3.0 or i = 9.0, respec- 
tively. For higher values of i, the approximation diverges 
quickly, just as the discrete gyroaverage method diverges 
from the analytic result for increasing k±pth- Based on 
this, care must be taken in selecting Ng for large values of 
k±pth- In this paper we restrict attention to k±pth< 3.0 
and use an adaptive Ng technique based on Ref. [29J . 



IV. DIAGNOSTICS FOR NON-DIFFUSIVE 
TRANSPORT 



called sub-diffusive. If 1 < 7 < 2, the spreading is faster 
than diffusion and transport is super-diffusive. A similar 
classification applies for sub-advection (0 < % < 1) and 
super-advection (1 < x < 2). In the presence of anoma- 
lous scaling, the introduction of an effective transport 
velocity or an effective diffusivity is meaningless since 
these transport coefficients are either zero (in the sub- 
advection/sub-diffusion case) or infinite (in the super- 
advection/super-diffusion case). The diagnostics based 
on the statistical moments are straightforward to imple- 
ment. The key is to look for a scaling region in a log-log 
plot of the moments as functions of time, after transients 
have passed. However, as with the data analyzed below, 
it is possible for the moments to follow different scaling 
regimes for different time intervals. 



In this section we review several Lagrangian diagnos- 
tics for transport study. After defining each diagnostic, 
we recall expected behavior for both diffusive and non- 
diffusive transport. These diagnostics have been success- 
fully used in transport experiments, models, and sim- 
ulations in both fluids and plasmas. For examples see 
Refs. [1^, [ia, [231 and references therein. To simplify the 
discussion we limit attention to one-dimensional trans- 
port, i.e. transport in the poloidal-like direction y. In 
the specific transport problem considered in this paper, 
y is in the direction of the propagation of the zonal flow 
and the drift waves, and is orthogonal to both the den- 
sity gradient and the magnetic field. Generalization of 
the diagnostics to higher dimensions is straightforward. 



A. Statistical moments of particle displacements 



The basic particle data consists of the ensemble 
{yi{t)}, with i = 1,2, ... Np, containing the time evo- 
lution of the y-coordinate of the Np test particles in the 
simulation. From here we define the ensemble of parti- 
cle displacements, {Syi{t)}, where Syi{t) — yi{t) — yi(0). 
The statistical moments of the particle displacements 
provide one of the simplest and most natural charac- 
terizations of Lagrangian transport. Of particular in- 
terest are the mean M{t) = {Sy) and the variance 
cr^(t) = {[Sy — {Sy)] ) where ( ) denotes ensemble aver- 
age. In the case of diffusive transport {e.g., a Brownian 
random walk), the moments exhibit asymptotic linear 
scaling in time, which allows the definition of an effec- 
tive transport velocity (pinch) V^ff and an effective dif- 
fusivity Deff according to Ve// = limj^cxj M{t)/t and 

However, in the case of nondiffusivc transport, the mo- 
ments display anomalous scaling of the form 



M r^t^. 



t^ . 



(15) 



with X 7^ 1 E^nd 7^1. If0<7<l the spread- 
ing is slower than in the diffusive case and transport is 



B. Particle displacement PDFs: spatial scaling 

The probability distribution function (PDF) of particle 
displacements, P{Sy,t\Sy',t'), contains all of the statisti- 
cal information from displacements beyond the first and 
second moments. By definition, P {Sy , t\Sy' , t' ^ t) ~ 
S{y). Numerically, P is constructed from the normalized 
histogram of particle positions at a given time. Formally, 
P{Sy, t\Sy', t') corresponds to the Green's function deter- 
mining the distribution of the test particles in terms of 
the initial particle distribution. For a Brownian random 
walk, the central limit theorem implies that P asymptot- 
ically approaches a Gaussian distribution, Pq, that sat- 
isfies diffusive scaling, Pq = t^^l'^G{Y jt^l"^), where G is 
a Gaussian and Y ~ Sy — {Sy). However, a non-diffusive 
propagator can exhibit the more general (anomalous) 
self-similar scaling 



P = t-^/^LiY/t'>^^) 



(16) 



where < 7 < 2 and L is a non-Gaussian function. Note 
that, by construction, the propagator has zero mean, 
and the scaling exponent 7 in Eq. (J16p is the same as 
the exponent in Eq. ([TS]). From Eq. (fTB|) it follows that 
P{Y, t) = X^''^P{X'/'^Y, \t) where A is a real number. 
Therefore, if the propagator is self-similar, P is invariant 
with respect to the space-time renormalization transfor- 
mation {Y, t) -^ {X^^^Y, Xt), up to a scale factor. 

Equation (J16p provides a useful diagnostic to reveal 
non-diffusive transport and, in particular, the existence 
of anomalous self- similar scaling. This diagnostic is im- 
plemented by plotting the propagator at different times 
in rescaled coordinates, i.e. t^^^P versus Y jVl"^ . With 
self-similar non-diffusive transport, the plots at different 
times rescale and collapse into a single function L. One 
of the most important departures from Gaussianity is al- 
gebraic decaying, "fat" tails in the propagator for large 
Sy at fixed t. 



P ^Sy 



-C 



(17) 



When this behavior is found, the value of the scaling 



exponent C is a useful diagnostic that characterizes the 
interniittency of the transport process. 



C. Trapping and flight probability distribution 
functions 

Diffusive transport can be interpreted as a coarse- 
grained (macroscopic) description of a fine-grained (mi- 
croscopic) Brownian random walk. In a similar way, non- 
diffusive transport can sometimes be viewed as the result 
of a non-Brownian random walk with a non-Gaussian 
and/or non-Markovian [30| underlying stochastic pro- 
cess. Trapping and flight probability distribution func- 
tions are two useful diagnostics for the characterization 
of non-Brownian random walks. Given a particle trajec- 
tory, Hiit), a trapping event is defined a portion of the 
trajectory during which the particle stays on an eddy. 
Flight events are portions that are not trapping events. 
Thus, each particle orbit in the ensemble of initial condi- 
tions may be decomposed as a sequence of trapping and 
flight events. 

Numerically, the events are detected by tracking re- 
versals in the Lagrangian acceleration of particles. From 
the histograms of trapping and flight events one may con- 
struct the probability distribution functions of trapping 
events, tp{t), and flight events, X{y). Indications of non- 
diffusive transport can be explored by studying the de- 
partures of X{y) and ijj{t) from the Gaussian and expo- 
nential dependencies characteristic of Brownian random 
walks. Of particular interest is the presence of asymp- 
totic algebraic scaling of the form, 

t/jr^r" , A - y~^ . (18) 

When /.J < 1 the mean waiting time, Jtipdt, is infinite 
and no characteristic temporal scale exists. In the Levy 
flight regime ^ < 3, and therefore the second moment, 
j if'Xdy, diverges and no characteristic spatial scale ex- 
ists. The PDFs of flight and trapping events are in prin- 
ciple interesting because of their connection to the con- 
tinuous time random walk (CTRW) model, which, in the 
fluid continuum limit, can be described using fractional 
diffusion equations 0, [3l|, l34| . 



D. Lagrangian velocity autocorrelation function 

Further insights into non-diffusive transport can be 
gained by looking at the Lagrangian velocity autocor- 
relation function C{t) = {vy(T)vy(0)) where Vy is the 
Lagrangian velocity of a particle. The Green-Kubo rela- 
tion, da^/dt = 2/„ C{T)dT, relates the velocity autocor- 
relation function to the variance of displacements. When 
C decays fast enough so that the integral converges, this 
relation can be used to deflne an effective diffusivity ac- 
cording to Deff = L C{T)dT. However, when C has 
algebraic decay of the form 



C{t) 



(19) 



with K < 1, the integral diverges and the concept of effec- 
tive diffusivity loses meaning. For super-diffusive trans- 
port, a^ ^ f implies 7 = 2 — k. 



V. NUMERICAL RESULTS 

For the Lagrangian statistics we consider ensembles of 
A^ = 8 X 10^ test particles, and integrate the equations of 
motion, with and without FLR effects, up to t = 5.2x 10^. 
The zero Larmor radius results were obtained from the 
numerical integration of the guiding center equations in 
Eq. ^ with the Hamiltonian in Eq. ^ with Si = 0.5, 
£2 = 0.2, ci = 0.4, kj_i = 6.0, fc_L2 = 5.0, uj = 6.0. 
The same Hamiltonian and parameter values were used 
in the FLR (0 < k±pth < 3) calculations based on an Ng 
adaptive gyroaverage. 

The Poincare plots in Fig. [3] show the dependence of 
the degree of stochasticity on the value of k±p. Fig- 
ure [3Ka) corresponds to k±p = 0. The degree of stochas- 
ticity is relatively large and, consistent with the results 
reported in Refs. [l|, [1^1, the stochastic layer is strongly 
asymmetric. In particular, the region of stochasticity left 
of the unperturbed separatrix (shown with the bold line) 
is very small. As will be discussed below, this asym- 
metry manifests directly in the skewness of the tail of 
the test particle propagator, which decays strongly for 
Sy > due to the very low probability of having sticky- 
flight particles jumping in the y > direction. It may be 
interesting to compare pth to the thickness of the lower 
branch of the stochastic region, A^. For example, when 
k_Lpth = {1.2,2}, Pth/ A, = {0.44,1.8}. This trend is 
mainly due to the rapid shrinkage of the stochastic layer 
as a function of ptt- When k±pth ~ 3, the value of Ag 
is very difficult to determine because the stochastic layer 
has almost completely disappeared. 

In the FLR calculations the test particles have a 
Maxwellian distribution of Larmor radii characterized by 
a mean value, pth. Thus, depending on its specific value 
of p, each particle "sees" a different Hamiltonian, which 
in general will be stochastic to a lesser degree as p in- 
creases. Figures [3ljb)-(d) illustrate this with Poincare 
plots corresponding to (b) k±p = 1.2, (c) k±p — 2.0 and 
(d) fcj_p = 3.0. Each one of these Poincare sections was 
computed by assigning the same value of fc_L/9, to all the 
initial conditions. It is observed that the value of k^p has 
a direct non-trivial influence on the degree of stochas- 
ticity. In general, a Poincare plot corresponding to an 
ensemble of particles with a Maxwellian distribution of 
gyroradii will be a mixture of k\^p Poincare plots, as seen 
in Fig. 31 The crossings of curves in the Poincare plots 
indicates the presence of multiple Hamiltonian systems 
indexed by values oik^p. 

To compute the Lagrangian diagnostics of non- 
diffusive transport, we considered groups of particles lo- 
cated in the vicinity of a hyperbolic equilibrium point of 
the Hamiltonian. The resulting trajectories can be di- 
vided into three categories: (a) passing trajectories that 
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FIG. 3: Dependence of phase space topology and stochasticity 
on Larmor radius for the Hamiltonian model in Eq. ((9)1. The 
panels show Poincare maps for a ensemble of particles with 
gyroradius distribution of the form H = S{k±p — k±pth) with 
(a) k±pth = 0, (b) k±pth = 1-2, (c) k±pth. = 2.0 and (d) 
k±pth = 3.0. The bold, solid curve indicates the unperturbed 
separatrix for k±pth = 0. 



FIG. 5: Typical sticky-flight trajectory in the Hamiltonian 
transport model. This particle alternates in a seemingly un- 
predictable way between being trapped in E x B eddies and 
being transported following the zonal shear flow. Other types 
of orbits, not shown, correspond to trapped orbits that never 
leave the original eddy, or passing orbits that move following 
the zonal flows without being trapped. 




FIG. 4: Poincare plot for multiple gyroradii values from the 
Maxwellian distribution with k±pth ~ 0.6. Crossings of 
curves indicate the presence of multiple Hamiltonian systems, 
one for each value of p. 



follow the zonal flow and never enter an E x B eddy 
(vortex), (b) stagnant trajectories which never leave an 
eddy and (c) sticky-flight trajectories which, as shown in 
Fig. [SI alternate between the eddies and the zonal flow. 
Since the statistics of the passing and the stagnant tra- 
jectories are trivial, these particles will be ignored during 
the data analysis. 

Several techniques for isolating sticky-flight trajecto- 



ries can be devised. Our trajectory filter works by ex- 
amining all trajectories during their entire history, and 
discarding those that never encircle a vortex (passing) 
and those that do not move more than one vortex width 
from ther original positions (stagnant). We have also 
tested a filter in Fourier-velocity space that discards hor- 
izontal velocity time series without a broadband spec- 
trum. Depending on the threshold for defining "broad- 
band," the Fourier filter gives practically the same re- 
sults as the trajectory filter. Analysis of sticky-flights in 
more realistic velocity flclds would be served better by a 
Fourier-velocity filter. The proper threshold for defining 
a "broadband" spectrum can be found from asymptotic 
considerations. 

Figure [S] shows the effect of the trajectory filter on the 
histogram of Larmor radii. In the computation of the his- 
togram we show the number of particles, A'^, multiplied 
by the appropriate metric factor p. The solid line denotes 
the histogram considering all the particles in the ensem- 
ble, i.e. without the filter. As expected, this histogram 
corresponds to a sampling of the Mawellian distribution 
in Eq. pT7)) . It is observed that the filter tends to remove 
particles with large p, and, as expected, the number of 
particles removed decreases with ti , the time of filter ap- 
plication. Since i; = 5200 appears to give an asymptotic 
value for the number of sticky-flights, it is used as the 
filtering time for the following diagnostics. When scaling 
values are reported for t < 5200, the filter is still applied 
uniformly at i = 5200. The first column in Table |T] gives 
Us, the percentage of sticky-flights, for each tested value 
of k±pth when the filter is applied at t; = 5200. 



^no filter (100%) 
^f, = 260 (6%) 
-^fi = 1040 (38%) 
-~U = 5200 (58%) 
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FIG. 6: (Color online) Gyroradius histogram for k±pth = 1-2 
with sticky-flight filter applied at various times. The upper- 
most curve shows the unfiltered distribution obtained from 
the sampling of the 2-D Maxwellian distribution in Eq. pO|l . 
The other curves give the distribution at different times af- 
ter the filter (which keeps only the sticky-flight orbits) has 
been applied. The vertical line marks the maximum of the 
unfiltered distribution. 



TABLE I: Measures of sticky trajectories and non-diffusive 
transport for the Vy = tanh(a;) model with initial positions in 
a box centered on a hyperbolic fixed point. The percentage 
of sticky trajectories at f = 5200, Ha, is shown, along with 
the mean and variance time power law exponents, x s^nd 7 
respectively, at early and late time. "Early" refers to a fit 
for 104 < t < 1040 and "late" refers to 4700 < t < 5200. 
Accuracy for these fits is similar to that observed in Fig. [T] 
and equal to ±0.1. 
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A. Super-diffusive scaling 

Before presenting the chaotic transport results, it is 
instructive to go back to the simple parallel flow in 
Eq. [11] to explore the role of FLR effects on particle 
dispersion in the context of an integrablc flow for a 



ensemble of particles initially distributed according to 
P = S{x — XQ)S{y — Ho). If all the particles have the 
same Larmor radius, i.e. if H(p) = 5{p — pth), then as 
Eq. IA3I in Appendix |^ shows, P maintains its delta 
function shape and simply drifts with the effective veloc- 
ity Jo{kp)Uo, which in the limit of zero Larmor radius 
corresponds to the parallel flow velocity. In this case, 
FLR effects are irrelevant since they simply rescale the 
velocity. However, when the particles have different Lar- 
mor radii, as in the Maxwellian case of Eq. IA41 the 
effective velocity of each particle will be different and the 
initial delta function will spread in space as is evident in 
the particle propagators shown in Fig.IH In this case, the 
first and second moments arc M = Vefft and a^ = At^, 
where Veff and A are functions of k±pth given in Ap- 
pendix 1^ The key issue to observe is that the variance 
docs not exhibit diffusive scaling, and that a distribution 
of Larmor radii gives rise to a ballistic spreading of the 
particles. 

For transport in the nonintegrable flow with the zonal 
flow and drift waves. Fig. [7] shows the mean, M{t), and 
variance, cr'^{t), for k±pth = and kj_pth = 0.6. A sum- 
mary of the values of the scaling exponents x ^^'^ 7 for 
all the values of k±pth studied is presented in Table [J To 
a good approximation, the mean exhibits linear scaling, 
i.e. X ~ 1 in Eq. (|15[) . indicative of regular advection, 
for all values of k±pth. The variance consistently shows 
clear evidence of super-diffusive transport, i.e. 7 > 1 
in Eq. psp . In the zero Larmor radius case, two scal- 
ing regimes are observed. Up to t « 10'^, which corre- 
sponds to the simulations in Ref. [l5[, the power law fit- 
ting in Fig.[7Ub) indicates an almost ballistic scaling with 
7 = 1.9. However, at a later time there is a transition 
to 7 = 1.6. As Table U shows, FLR effects seem to elim- 
inate the distinction between early and late regimes. In 
particular, according to Fig. [7]^d) where k±pth = 0.6, the 
scaling 7 = 1.6 holds throughout the integration time. 
As a general trend, it is observed that the exponent 7 
decreases with increasing k±pth beyond 0.1. Statistics 
for sticky-flights become poor for k±pth = 3 because the 
degree of stochasticity [see Fig. ^d)] becomes small. 



B. Asymmetric, non-Gaussian PDF of particle 
displacements 

Motivated by the presence of two different scaling reg- 
imens in the variance, we study the PDF of particle dis- 
placements at intermediate and large times. Figure [5] 
shows the PDFs at intermediate times, with [HJa) cor- 
responding to k±pth ~ and [UJb) corresponding to 
k±pth = 1.2. The solid lines denote the PDFs of the fil- 
tered data, (i.e. including only sticky-flight orbits) and 
the dashed line denotes the PDFs of the unfiltered data. 
The spikes for large Sy in the unfiltered distributions re- 
sult from the contribution of passing orbits that the fil- 
ter effectively removes. The filtered PDFs are clearly 
non-Gaussian with strong skcwness in the negative Sy 




FIG. 7: (Color online) Time evolution of statistical moments 
of particle displacements. Panels (a) and (b) correspond to 
k±pth = and panels (c) and (d) correspond to k±pth = 0.6. 
Plots (a) and (c) give the absolute value of the first moment 
M, and plots (b) and (d) show the second moment. The 
dashed lines in panels (a) and (c) have slopes corresponding 
to )(; = 1.1(0.9) and x = 1-0 indicative of normal advection 
scaling, i.e. |M| ~ f-^ with x ~ 1- The variance shows super- 
diffusive scaling i.e. a^ ~ f with 7 7^ 1. However, in the 
k±pth = case, a sharp transition is observed in the anoma- 
lous diffusion exponent. The dashed lines in panels (b) have 
slopes corresponding to 7 = 1.9 and 7 = 1.6. The dashed line 
in panel (d) has a slope corresponding 7 = 1.9 indicating a 
uniform scaling of the variance for k±pth = 0.6. 



direction. The strong left-right asymmetry of the PDFs 
results from the asymmetry of the stochastic layer. 

In particular, as the Poincare plots in Fig. [3] show, 
the stochastic layer is thicker on the right side of the 
vortex. This asymmetry depends on the value of the 
perturbation frequency u in Eq. ([9]). In fact, as discussed 
in Ref. [l^, the relative thickness of the stochastic layers, 
and therefore the symmetry of tracer transport, can be 
controlled by changing uj. As the insets in Fig. 7 show, 
both PDFs decay algebraically as in Eq. ([17]). However, a 
strong dependence of the decay exponent on the value of 
the Larmor radius is observed. For k±pth = 0, ^ ~ 1.95, 
and for k±pth = 1-2, C ~ 2.9. As Table |T] indicates, the 
value of the decay exponent C. increases monotonically 
with k±pth- 

The particle displacement PDFs at longer times are 
shown in Fig. [HI As before, the solid lines denote the 
filtered distribution and the dashed lines the unfiltered 
distribution. A critical dependence on the Larmor radius 
is observed. For k±pth — the PDF transitions to an 
exponential decaying distribution, whereas for k±pth = 
0.6 the PDF maintains its algebraic decay with the same 
exponent as the one observed at short times, C « 2.9. 
The robustness of the algebraic decay in the finite Larmor 
radius case might be attributed to the persistence of large 




FIG. 8; (Color online) Probability distribution function of 
particle displacements at intermediate times, t = 1040. Panel 
(a) corresponds to k±pth = and panel (b) corresponds to 
k±pth ~ 1.2. The insets in both figures show evidence of 
algebraic decaying tails, P ~ (5j/~'' with (" — 1.95 for k±pth = 
and ( = 2.9 for k±pth ~ 1.2. In both plots, the solid line 
denotes the PDF of sticky-fiights {i.e., excluding the passing 
and trapped orbits), and the dashed line denotes the PDF 
computed using all the orbits. 



particle displacements which, due to the presence of the 
strong zonal flows, are enhanced by the gyroaverage. One 
should note that a Levy process requires C < 3, which 
means that the increase of k±pth moves the process away 
from the Levy type. 

The transition from algebraic to exponential decay in 
the zero Larmor radius case is likely due to the presence 
of truncated Levy flights. Exact Levy flights produce 
long particle displacements that result in slowly decay- 
ing algebraic tails at all times. However, non-ideal ef- 
fects such as particle decorrelation might preclude the 
existence of arbitrarily long displacements, resulting in a 
faster than algebraic decay of the tails at long times. See, 
for example, Refs. [33, [SJ, [3a] for more details on trun- 
cated Levy processes. One obvious reason for a truncated 
Levy process in the present system is the finite veloc- 
ity requirement, which precludes the existence of infinite 
jumps. 



C. Levy flights and algebraic trapping PDFs 

Figure [TO] shows the trapping time and flight length 
PDFs for k±pth = in (a) and (c), and for k±pth = L2 
in (b) and (d). In both cases, the trapping PDF clearly 
decays algebraically as in Eq. ^TE\i . with i^ = 1.8 for 
k±pth = 0, and ly = 2.0 for k±pth ~ 1.2. Figures [TOTc) 
and ^d) show the PDFs of flight lengths. Note that, 
because transport in this case is asymmetric, there are 
actually two flight PDFs, one corresponding to positive 
flights (with dashed fit line) and another corresponding 
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FIG. 9: (Color online) Probability distribution function of 
particle displacements at large times, t = 5200. Panel (a) cor- 
responds to k_i_pth = and panel (b) corresponds to k±pth ~ 
1.2. In case (a) the PDF decays exponentially, P ~ e""^ " 
with A ~ 0.002. On the other hand, for k±pth = 1.2, 
the inset shows evidence of algebraic decay, P ~ Sy~'' with 
^ = 2.9. In both plots, the solid line denotes the PDF of 
sticky-flights {i.e., excluding the passing and trapped orbits), 
and the dashed line denotes the PDF computed using all the 
orbits. 



to negative flights (solid fit line). The PDF of negative 
flights decays as a power law with ^ = 1.8 for k±pth = 0, 
and pL = 2.7 for k±pth = 1-2. Since /i < 3 in both cases, 
these flights correspond to Levy flights. However, the de- 
cay of the curve for positive flights is much steeper with 
/^ > 3 regardless of the value of k±pth, which implies that 
positive displacements are not Levy flights. The tails of 
the trapping and flight PDFs transition to exponential 
decay at SyfUght ~ —1000 and ttrapt ~ 2000. As dis- 
cussed before, this transition is indicative of the possible 
presence of truncated Levy flights. 



D. Algebraic decay of Lagrangian velocity 
autocorrelation function 



Figure [TT] shows the Lagrangian velocity autocorrela- 
tion function for the sticky-flights with k±pth = in 
Fig. [nja) and with k^pth == 1-2 in Fig. \Tllh). Both 
curves follow algebraic decay of the form C(t) ^ r"". 
When k±pth = 0, k = 0.2 and when k±pth = 1.2, k = 0.3. 
Both values are consistent with the Green-Kubo relation 
between the decay of the velocity correlation and the 
scaling of the variance according to which k — 2 — 7. 
The frequency of small scale oscillations observed in the 
correlation seems to increase when k±pth changes from 
0^ 1.2. 



FIG. 10: (Color online) Probability distribution functions of 
particle trapping events and particle flight events for k±pth = 
and k±pth ~ 1.2. The trapping PDFs are shown in (a) and 
(b), and the flight PDFs are shown in (c) and (d). Panels 
(a) and (c) correspond to k±pth ~ 0, and panels (b) and 
(d) correspond to k±pth ~ 1-2. The solid straight lines in 
(a) and (c) indicate that the trapping PDFs show algebraic 



decay, P >- 
for k±pth 



with v 



for k\ 



0, and ty « 2.0 



trap! vviuii 1^ ~ i.u iui i^^pth 

1.2. The negative flights PDF shown fit with 



solid lines also exhibit algebraic decay of the form P ~ t~fuqht 
with p « 1.8 for the case k±pth = 0, and p Ri 2.7 for the 
case k±pth = 1.2. The PDFs of positive flights, shown fit 
with dashed lines, show a faster exponential-type decay with 
p « 3.0 in both cases. 
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FIG. 11; Lagrangian velocity autocorrelation function for 
sticky-flight trajectories. Panel (a) corresponds to k±pth = 
and panel (b) corresponds to k±pth = 1.2. The curves with 
dots are the numerical results, and the solid line curves are 
algebraic fits of the form C ~ r"" with k = 0.2 in (a) and 
K = 0.3 in (b). 
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VI. SELF-SIMILAR ANOMALOUS SCALING 
AND FRACTIONAL DIFFUSION MODELING 

An important goal of transport modeling is to con- 
struct effective transport equations that describe the 
"macroscopic" coarse grained dynamics when given in- 
formation at the "microscopic" kinetic level. When 
the microscopic dynamics involves Gaussian, Markovian 
stochastic processes {e.g., a Brownian random walk) the 
macroscopic dynamics can be modeled using diffusion 
type equations. This is the basic idea behind the use of 
diffusive models to describe coUisional transport. How- 
ever, in recent years it has been shown that the standard 
diffusion picture can fail when non- Gaussian and/or non- 
Mar kovian statistics are present. 

In particular, experimental, numerical and analytical 
transport studies in fluids and plasmas {e.g. Refs. [ll|, 
E, m M, m, S M, m, S \^ and references therein) 
have shown that underlying stochastic processes govern- 
ing particle transport in flows with coherent structures, 
like zonal flows and eddies, typically involve anomalously 
large particle displacements induced by the zonal flows 
and/or anomalous particle trapping in eddies. The pres- 
ence of large particle displacements can invalidate the 
Gaussianity of displacement distributions. Particle trap- 
ping can introduce waiting time effects that invalidate the 
Markovian assumption because of memory effects. The 
statistics of particle transport discussed in the previous 
section shows clear evidence of these type of phenomena. 
This section presents an effective macroscopic model that 
describes quantitatively the spatio-temporal evolution of 
the PDF of particle displacements. 

An important piece of information needed for con- 
structing an effective transport model is shown in 
Fig. [T21 Figures [T^ a)-(c') show the temporal evolution 
of the PDF of particle displacements for different values 
of k±pth- As discussed before, the PDF develops a strong 
"fat" tail to the left and, by conservation of probability, 
the peak of the distribution goes down. Figures [TWdl-ff) 
show the same data plotted using rescaled variables as in 
Eq. P^ . In the horizontal axis, 77 = Sy/f^^^, and in the 
vertical axis, P has been multiplied by the factor t"^'^, 
where 7 is the anomalous diffusion exponent in Eq. psp . 
From this it follows that the PDF at a time Xt is re- 
lated to the PDF at time t by the scaling transformation 
P{Sy,Xt) = A-T/2p(y/A^/2,i). The fact that, for the 
problem of interest here, 7 7^ 1, rules out the possibility 
of constructing a transport model based on the diffusion 
equation with an effective diffusivity because the solution 
of the diffusion equation scales as P = t~^^'^L{dy/t^^^). 

A natural way to built transport models that display 
self-similar anomalous scaling is to use fractional diffu- 
sion equations of the from 
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FIG. 12: (Color online) Self-similar scaling of probability 
distribution function of particle displacements (PDF). The 
curves denote the PDFs at i = 1040, t = 1560, and t = 2080, 
with later times showing more spreading in the PDF. Panels 
(a), (b) and (c) show the PDFs corresponding to k±pth = 0, 
k±pth = 0.6 and k±pth = 1.2, respectively. Panels (d), (e) 
and (f) show the collapse of the corresponding PDFs when 
plotted as functions of the similarity variable 77 = Sy/t"'''^ 
and rescaled with the factor t^". 



left and right fractional derivatives. These non-local op- 
erators are a natural generalization of the regular dif- 
ferential operator, 9" of integer order n. For example, 

Fourier transforms of the fractional operator, Tlf] = f = 
J e'-'^^fdy, satisfy 

^ [-ooD'^p] = (-^fc)" p ' ^ [yD^p] = a^r P , 

(21) 
for non-integer values of a. In a similar way, the operator 
on the left hand side of Eq. (PU)) is a natural extension 
of the regular time derivative, dtf , in the sense that its 



Laplace transform, £[/] = f = J e '"'^fdt, satisfies 



L ^D^P 



s^P- 



s^-^P{t^O). 



(22) 







D'^P^Xf[l-ooD'^+ryD^]P, 



(20) 



for < /? < 1. As expected, Eq. (^0]) reduces 
to the standard diffusion equation when a = 2 and 
(3 = 1. Further formal details on fractional deriva- 
tives, including their representation in the y and t do- 
mains in terms of non-local operators can be found in 
Refs. [4l|, I44I . For a discussion on the use of these opera- 
tors to model non-diffusive transport in plasmas, see for 
example Refs. [T^ |T^ and references therein. 

To explore the self-similarity properties of the frac- 
tional diffusion model we use Eqs. (PT|) - ([^^ and write 

the Fourier-Laplace transform, Q, of the Green's func- 
tion, g, of Eq. (HOI) as 



where/ — — sec(a7r/2)(l — 0)/2, and r = — sec(Q;7r/2)(l-|- 
0)/2. The operators -00-Dy and yD"^ are called the 






A = Xf[n-^kr+r{ikr] 



(23) 
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where a ^ 1 and G{y,t = 0) = S{y). It follows di- 
rectly from Eq. ^ that g{k, s/X) = Xg{X^^°'k, s) which 
in y-t space implies the self-similar scaling g{y, Xt) — 
^-/V«g(^^-/9/"y^^) of the fractional diffusion propagator 
Eq. (|20p . Therefore, the fractional equation will exhibit 
the same self-similar scaling as the numerically obtained 
PDF provided the fractional orders of the spatial and 
temporal derivatives satisfy 



7 = 2/3/a . 



(24) 



According to Table 1, to a good approximation, 7 « 2 
in the intermediate asymptotic regime. Based on this ob- 
servation, and following Eq. (j24p . we will assume a = /? in 
the fractional diffusion model. This special case, known 
as neutral fractional diffusion, has a Green's function that 
can fortunately be expressed in closed form using elemen- 
tary functions, as shown in Ref. [isj : 



g{v;a,§) = - 



sin TT{a~9)/2 rj' 



^ 1 + 2?7" cos 7r(a - e)/2 + if"" 

for 77 > , (25) 

where 77 = Sy/fi"^ is the similarity variable and \6\ < 
min{a, 2 — a}. The solution for 77 < is obtained using 
the relation g{—ri;a,9) = g{T];a,—0). The parameter 
is related to the asymmetry parameter 9 introduced be- 
fore in the definition of the weighting factors I and r ac- 
cording to = tan(7r0/2)/tan(7ra/2). Given the Green's 
function, the solution of the fractional diffusion equation 
for an initial condition Po{Sy) = P{Sy,0) is 



Pi6y, t) = / P^{5y') g{5y - 6y', t)dSy' 



(26) 



For the initial condition we assume a localized distri- 
bution of the form P) = 1/A for \dy\ < A/2 and Pq = 
elsewhere (see Ref. [I^ ) . The use of this initial condition 
is necessary to account for the presence of transients in 
the evolution of the PDF not reproduced by the fractional 
diffusion equation, which describes the intermediate time 
regime. Figure [13] shows the comparison of the solution 
of the fractional diffusion equation in Eq. (PO)) according 
to Eqs. (P5|) and (PS)) and the numerically obtained PDF 
obtained from the histograms of particle displacements 
at i = 936 for k±pth — in Fig. ll3r a) and k±pth ~ 0.6 in 
Fig. [TWb) . For the fractional diffusion model parameters 
we used a ~ P — 0.80 and 6 = 0.79 in the k±pth = case, 
and a — /3 — 0.85 and 6 = 0.84 in the k±pth = 0.6 case. 
In both cases, we used A = 60, which is small compared 
to the maximum range of the PDF, Sy ^ —800. 



VII. SUMMARY AND CONCLUSIONS 

In this paper we presented a numerical study of FLR 
effects on non-diffusive transport of test particles in a flow 
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FIG. 13: (Color online) Comparison between the numerically 
determined PDF of particle displacements (square markers) 
and the solution of the effective fractional diffusion model in 
Eq. (|25|) (solid lines) with a localized initial condition. In panel 
(a), which corresponds to k±pth = 0, the fractional diffusion 
parameters are a — 0.8, 9 = 0.79, yl = 60 and Xf = 0.15. 
For the case k±pth ~ 0.6, shown in panel (b), a — ft — .85, 
9 = 0.84, A = 60 and X = 0.12. 



dominated by a strong zonal shear flow and large scale 
E X B eddies. We modeled the flow using a Hamiltonian 
dynamical system consisting of a linear superposition of a 
strong zonal shear flow and eigenmodes of the HM equa- 
tion. For the parameter values considered, the Hamil- 
tonian causes chaotic transport. Test particles alternate 
stochastically between being trapped in the vortices and 
being transported by the zonal flow. To expose the non- 
diffusive properties of the system we used Lagrangian 
statistical diagnostics including: (i) the scaling in time 
of statistical moments; (ii) the PDFs of particle displace- 
ments, (iii) trapping events and (iv) flight events; and 
(v) the decay of the Lagrangian velocity autocorrelation 
function. 

Finite Larmor radius effects were incorporated in the 
particle calculations by substituting the value of the 
E X B velocity at the location of the guiding center by 
its value averaged over a ring of radius p, where p is the 
Larmor radius. The ring average was computed using 
a discrete approximation. The numerical method was 
benchmarked using an analytical solution for a parallel 
zonal flow with no waves. We found that for k±p < 3 an 
8-point average gives accurate results, but higher order 
approximations must be used for for k±p > 3. Con- 
trary to previous works where all the particles were as- 
sumed to have the same value of p, here we considered 
a more realistic Maxwellian distribution of Larmor radii. 
Poincare plots revealed that the Larmor radius has a di- 
rect nontrivial effect on the topology of the flow and the 
degree of chaos of test particles. In particular, it was ob- 
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served that the amount of chaos, measured by the width 
of the stochastic region, is significantly reduced as fc^pj^ 
increases from to 3. A distribution of Larmor radii 
can also have a direct effect on the dispersion of parti- 
cles. In particular, we have shown that, even in the case 
of a completely integrable flow, particles exhibit ballis- 
tic spreading, cr^ ~ i^, when they have different Larmor 
radii. 

For the Lagrangian statistics we limited attention to 
sticky-flight orbits and ignored trapped and passing or- 
bits. The rationale for this fllter is that the trivial dynam- 
ics of passing and trapped particles give rise to outliers 
that artificially bias the statistics. The first moment, to 
a good approximation, has normal advective scaling, i.e. 
M ^ t^^ with X ~ 1, and the second moment has super- 
diffusive scaling, i.e. a^ r^ H, with 7 > 1. For k^p = 0, 
a sharp transition was observed in the scaling exponent, 
from 7 = 1.9 at intermediate times to 7 = 1.6 at larger 
times. Similar transitions in the value of 7 have been 
also found in other systems including temporally irregu- 
lar channel flows [36|, time dependent, three dimensional 
flows [4J], and two-dimensional vortex flows [lg|. For 
specific experimental instances, early time behavior will 
be more important than late time behavior if the do- 
main crossing time is small enough. We have found that 
FLR effects seem to eliminate the distinction between 
early and late time. For the range of k±pth considered, 
7 w 1.8 ± 0.1. We refer to this regime as super-diffusive 
ballistic transport since the variance approaches ballistic 
scaling (7 = 2) but the PDF of displacements retains a 
super- diffusive a ppe arance. Complementary results were 
obtained in Ref. [l2l | for nonlinear HM simulations. 

We also observed that the Lagrangian velocity auto- 
correlation function decays algebraically, C ^ t^'' where, 
in reasonable agreement with the Grecn-Kubo scaling, 
C = 2 — 7. The trapping and flight distributions show 
algebraic decay. The trapping time exponent, v, remains 
the same when k±pth changes. The PDFs of negative 
flights qualify as truncated Levy distributions but posi- 
tive flights are definitively not Levy. The negative flight 
exponent for k±pth = 1.2 is larger than expected in the 
context of a CTRW. 

At intermediate times, consistent with Refs. [i|,[ia|, the 
PDF of particle displacements in the zero Larmor radius 
case is an asymmetric non-Gaussian distribution with 
an algebraic decaying leftward tail. However, for larger 
times, the tail of the PDF transitions from algebraic to 
exponential decay. This algebraic-exponential transition 
in the PDF is likely to be related to the presence of trun- 
cated Levy flights, which, as discussed in Ref. [35|, might 
result from particle decorrelation or the finite size of pos- 
sible displacements. The robustness of the algebraic de- 
cay in the finite Larmor radius case might be attributed 
to the persistence of large particle displacements which, 
due to the presence of the strong zonal flows, are en- 
hanced by the gyroaverage. We have also shown that 
the PDF of particle displacements has self-similar scal- 
ing behavior for < k±pth < 3 and k±pth 7^ 0. Most 



importantly, we have shown that these distributions cor- 
respond to solutions of the neutral (a = (3) asymmetric 
fractional diffusion equation. 

Future work will apply the ideas and tools developed 
here to turbulent flows to more realistic plasma tur- 
bulence models. In particular, we will examine self- 
consistent particle transport parallel to a density gradi- 
ent in a gyrokinetic particle-in-cell simulation. Transport 
properties of tracers and self-consistent particles should 
be compared. 
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APPENDIX A: GYRO-AVERAGED PARTICLE 
PROPAGATOR IN A PARALLEL FLOW 

The gyroaverage equations of motion for test particles 
in the parallel zonal flow of Eq. [TT] are 



dx 



— = -0o^i (sin(fc_La;))8 = 

— (/^o^-L'^lfciP) sin(fcj^x) . (Al) 



A straightforward integration assuming an intial condi- 
tion (a;o,2/o) gives 

x = Xf), y = yo-UoJo{k±p)t, (A2) 

where Uq — (j)ok±sui{k±xo). From here it follows that 
the two-dimensional propagator is 

rir,t\r\t';p)^Six-x')6[y-y' + .Mk^p)Uot] . (A3) 

Integrating over x and assuming a Maxwellian distribu- 
tion of gyroradii gives the one-dimensional propagator in 

y, 

P{y,t\y',t';p) = 

) 

S[y-y' + Jo{k^p)Uot]pe-p"'p"<^ dp . (A4) 



PlhJo 



Integrating over p using basic properties of the delta func- 
tion gives Eq. (T^l From Eq. (jA4p it follows that the n-th 
moment of the gyrocenter displacement 5y ~ y — y' scales 
like t" according to 



{{SyD^iUotr 4\k^p)H{p)dp. (A5) 

Jo 
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where H{p) is the gyroradii distribution function. For 
n = 1 and n = 2 we recover the moments in Sec. FVTA) 
with 



v; 



eff 



t/ne-'^->-/4 



A 



U^e 



-felpL/2 



[lo (fclpL/2) 



(A6) 



in the case when H is Maxwelhan, where /q is the modi- 
fied Besscl function of zero-order. It is interesting to note 
that A has a maximum for k±pth ~ 2.5. 



[1] D. del Castillo Negrete, Phys. Plasmas 7, 1702 (2000). 
[2] A. Hasegawa and K. Mima, Physics of Fluids 21, 87 [26 
(1978). [27 

[3] W. W. Lee, Journal of Computational Physics 72, 243 

(1987). [28 

[4] R. Metzler and J. Klafter, Physics Reports 339, 1 (2000). 
[5] J.-P. Bouchaud and A. Georges, Physics Reports 195, [29 

127 (1990). 
[6] R. Balescu, Aspects of Anomalous Transport in Plasmas [30 

(lOP Puslishing, Bristol, 2005). 
[7] D. del Castillo-Negrete, in Turbulent transport in fusion [31 
plasmas: First ITER International Summer School , 
edited by S. Benkadda (AIP, College Park, 2008). [32 

[8] W. Horton, Plasma Physics 23, 1107 (1981). 
[9] G. Manfredi and R. Dendy, Phys. Rev. Lett. 76, 4360 
(1996). 

[10] G. Manfredi and R. Dendy, Phys. Plasmas 4, 628 (1997). [33 

[11] S. Benkadda, P. Gabbai, and G. M. Zaslavsky, Phys. 

Plasmas 4, 2864 (1997). [34 

[12] S. Annibaldi, G. Manfredi, R. Dendy, and L. Drury, [35 

Plasma Phys. Control. Fusion 42, L13 (2000). 
[13] S. V. Annibaldi, G. Manfredi, and R. O. Dendy Phys. [36 

Plasmas 9, 791 (2002). 
[14] W. Horton and A. Hasegawa, Chaos 4, 227 (1994). [37 

[15] D. del Castillo Negrete, Phys. Fluids 10, 576 (1998). 
[16] S. Kovalyov, Chaos 10, 153 (2000). [38 

[17] A. Chechkin, V. Gonchar, and M. Szydlowski, Phys. 

Plasmas 9, 78 (2002). [39 

[18] D. del Castillo Negrete, Phys. Plasmas 11, 3854 (2004). [40 

[19] D. del Castillo Negrete, Phys. Plasmas 13, 082308 (2006). 
[20] L. Garcia and B. A. Carreras, Physics of Plasmas 13, [41 

022310 (2006). 
[21] D. del Castillo-Negrete et al, Nuclear Fusion 48, 075009 [42 

(13pp) (2008). 
[22] I. Calvo et al, Physics of Plasmas 15, 042302 (2008). 
[23] F. Jenko and W. Dorland, Phys. Rev. Lett. 89, 225001 [43 

(2002). 
[24] I. Furno et al, Physical Review Letters 100, 055004 [44 

(2008). 
[25] T. Hauff, F. Jenko, and S. Eule, Physics of Plasmas 14, 



102316 (2007). 

H. Aref, Journal of Fluid Mechanics 143, 1 (1984). 
T. Solomon, E. R. Weeks, and H. L. Swinney, Physical 
Review Letters 71, 3975 (1993). 

J. M. Finn and D. del Castillo Negrete, Chaos 11, 816 
(2001). 

A. Mishchenko, A. Konies, and R. Hatzky, Phys. Plasmas 
12, 062305 (2005). 

M. Shlesinger, G. M. Zaslavsky, and J. Klafter, Nature 
(London) 31, 363 (1993). 

E. W. MontroU and G. H. Weiss, Journal of Mathemati- 
cal Physics 6, 167 (1965). 

E. W. MontroU and M. F. Shlesinger, in Nonequilibrium 
Phenomena II. From Stochastics to Hydrodynamics , 
edited by J. L. Lebowitz and E. W. MontroU (Elsevier, 
Amsterdam, 1984). 

R. N. Mantegna and H. E. Stanley, Physical Review Let- 
ters 73, 2946 (1994). 

I. Koponen, Phys. Rev. E 52, 1197 (1995). 
A. Cartea and D. del Castillo-Negrete, Physical Review 
E 76, 041105 (2007). 

S. C. Venkataramani, T. M. Antonsen, Jr., and E. Ott, 
Physical Review Letters 78, 3864 (1997). 

E. R. Weeks and H. L. Swinney Phys. Rev. E 57, 4915 
(1998). 

F. Dupont, R. I. McLachlan, and V. Zeitlin, Phys. Fluids 
10, 3185 (1998). 

T. Benzekri et al, Phys. Rev. Lett. 96, 124503 (2006). 

S. V. Prants, M. V. Budyansky, and M. Y. Uleysky, 

Chaos 16, 033117 (2006). 

I. Podlubny, Fractional Differential Equations (Academic 

Press, San Diego, 1999). 

S. G. Samko, A. A. Kilbas, and O. I. Marichev, Fractional 

Integrals and Derivatives: Theory and Applications 

(Taylor and Francis Books Ltd, London, 1993). 

F. Mainardi, Y. Luchko, and G. Pagnini, Fract. Calc. 

App. Analysis 4, 153 (2001). 

M. A. Fogleman, M. J. Fawcett, and T. H. Solomon, 

Physical Review E 63, 020101 (2001). 



